import ij.IJ;
import ij.ImagePlus;
import ij.gui.NewImage;
import ij.io.FileSaver;
import ij.plugin.PlugIn;
import ij.process.ColorProcessor;
import ij.process.ImageProcessor;


import java.awt.Color;


import java.awt.Font;

import java.awt.Polygon;
import java.io.PrintStream;

import java.util.ArrayList;
import java.util.Iterator;
import java.util.Locale;
import java.util.Random; 



import pluginUtils.Circle;
import pluginUtils.CircleArrangement;
import pluginUtils.Complexe; 
import pluginUtils.DoyleCircle;
import pluginUtils.Droite;
import pluginUtils.MobiusTransform;
import pluginUtils.CircleInversion;
import pluginUtils.PackCircle;
import pluginUtils.Plane;
import pluginUtils.Point2D;
import pluginUtils.Pos3D;
import pluginUtils.Sphere;
import pluginUtils.SphereInversion;
import pluginUtils.Vector2D;
import sliceformUtils.CoupleDAngles;
import sliceformUtils.DroiteSpace3D;
import sliceformUtils.Ligne;
import sliceformUtils.Mesh;
import sliceformUtils.STLReader;
import sliceformUtils.Saucisson;
import sliceformUtils.Segment;
import sliceformUtils.Triangle;


public class CyclideDupin extends ImagePlus implements PlugIn{
	{Locale.setDefault(Locale.US); }
	
	
	private double mu=2; // r0 
	private double c=1;  // r1
	private double a=4;  // rc
	private double b=Math.sqrt(a*a-c*c); 
	
	private int nbCercles=10; 
	private int nbDiv=10; // toujours pair ? 
	
	
	private ImageProcessor primeIP; 
	private int X=2000,Y=2000; // image size
	
	private PrintStream output;
	private double xShift=0; 
	private double zShift=0; 
	
	
	
	protected DroiteSpace3D[][] intersections; 
	
	private double calculX(double th, double al) {
		double den=a-c*Math.cos(th)*Math.cos(al); 
		return (mu*(c-a*Math.cos(th)*Math.cos(al))+b*b*Math.cos(th))/den; 
	}
	
	private double calculY(double th, double al) {
		double den=a-c*Math.cos(th)*Math.cos(al); 
		return b*Math.sin(th)*(a-mu*Math.cos(al))/den; 
	}

	private double calculZ(double th, double al) {
		double den=a-c*Math.cos(th)*Math.cos(al); 
		return b*Math.sin(al)*(c*Math.cos(th)-mu)/den; 
	}

	private Random generator=new Random(); // One always needs some randomness
	
	// Scaling  the image
	private double x0=0,x1=1,y0=0,y1=1;
	private double acoef,ccoef,vcoef,wcoef; 
	
	public void setCoefs(){
		acoef=X/(x1-x0); 
		vcoef=Y/(y1-y0); 
		ccoef=-acoef*x0; 
		wcoef=-vcoef*y0; 
	}
	 
	
	
	public void setBounds(double xl,double xu,double yl,double yu){
		x0=xl; 
		x1=xu;
		y0=yl;
		y1=yu; 
	}
	
	private void drawPoint(double x,double y,ImageProcessor ip){
		int intX=(int)(x*acoef+ccoef);
		int intY=(int)(y*vcoef+wcoef);
		ip.fillOval(intX-10,intY-10,20,20);
	}
	
	private void drawLine(double xd,double yd,double xf,double yf){
		int XD=(int)(xd*acoef+ccoef); 
		int YD=(int)(yd*vcoef+wcoef); 
		int XF=(int)(xf*acoef+ccoef); 
		int YF=(int)(yf*vcoef+wcoef);
		primeIP.drawLine(XD,YD,XF,YF); 
	}
	
	private void drawLine(double xd,double yd,double xf,double yf,ImageProcessor ip){
		int XD=(int)(xd*acoef+ccoef); 
		int YD=(int)(yd*vcoef+wcoef); 
		int XF=(int)(xf*acoef+ccoef); 
		int YF=(int)(yf*vcoef+wcoef);
		ip.drawLine(XD,YD,XF,YF); 
	}
	
	private void drawLine(Point2D a,Point2D b,ImageProcessor ip){
		drawLine(a.getX(),a.getY(),b.getX(),b.getY(),ip); 
	}
	
	private void drawCircle(Vector2D center,double rayon,ImageProcessor ip){
		double xCenter=center.getX(); 
		double yCenter=center.getY(); 
		int intX=(int)((xCenter-xShift)*acoef+ccoef);
		int intY=(int)((yCenter-zShift)*vcoef+wcoef);
		int intRx=(int)(rayon*acoef);
		int intRy=(int)(rayon*vcoef);
		ip.drawOval(intX-intRx,intY-intRy,intRx*2,intRy*2);
	}
	
	public ArrayList<Pos3D> intersectCadre(double ymin,double ymax,double zmin,double zmax,DroiteSpace3D Q){
		ArrayList<Pos3D> resu=new ArrayList<Pos3D>(); 
		Pos3D A=new Pos3D(0,ymin,zmin); 
		Pos3D B=new Pos3D(0,ymin,zmax); 
		Pos3D C=new Pos3D(0,ymax,zmin); 
		Pos3D D=new Pos3D(0,ymax,zmax); 
		if(coteVsDroite(A,B,Q)!=null) resu.add(coteVsDroite(A,B,Q)); 
		if(coteVsDroite(A,C,Q)!=null) resu.add(coteVsDroite(A,C,Q)); 
		if(coteVsDroite(B,D,Q)!=null) resu.add(coteVsDroite(B,D,Q)); 
		if(coteVsDroite(C,D,Q)!=null) resu.add(coteVsDroite(C,D,Q)); 
		return resu; 
	}
	
	/* Trouver l'intersection de la droite avec un des quatres cotes de l'image
	 * retourne null si parallele ou en dehors de l'image. 
	 * on ne s'occupe pas de la valeur de x (tout se passe dans le plan perpendiculaire
	 * a l'axe Ox)
	 */
	public Pos3D coteVsDroite(Pos3D A,Pos3D B,DroiteSpace3D D){
		double yv=D.getDirection().getY(); 
		double zv=D.getDirection().getZ(); 
		double ya=A.getY(); 
		double za=A.getZ(); 
		double yb=B.getY(); 
		double zb=B.getZ(); 
		double y0=D.getPoint().getY(); 
		double z0=D.getPoint().getZ(); 
		double det=-zv*(ya-yb)+yv*(za-zb); 
		if(Math.abs(det)<1e-6) return null; 
		double detAlpha=-zv*(y0-yb)+yv*(z0-zb);
		double alpha=detAlpha/det; 
		if ((alpha<0)||(alpha>1)) return null; 
		// on s'en fout de X, on va tracer une droite dans un plan
		double ny=alpha*ya+(1-alpha)*yb; 
		double nz=alpha*za+(1-alpha)*zb;
		return new Pos3D(0,ny,nz); 
	}
		
	public DroiteSpace3D inter(Plane p1,Plane p2){
		if((p1.getA()==0)&&(p2.getA()!=0)){
			Plane temp=p1; p1=p2; p2=temp;
		}
		if((p1.getA()==0)&&(p2.getA()==0)){ // droite de direction (1,0,0) (ou bien plan paralleles, pas de solution)
			if((p1.getB()==0)&&(p2.getB()==0)) return null; // plans paralleles
			if((p1.getB()==0)&&(p2.getB()!=0)){
				Plane temp=p1; p1=p2; p2=temp;
			}
			// Chercher un point (x,y,z)
			double bdm=p2.getB()*p1.getD()-p1.getB()*p2.getD(); 
			double cbm=p2.getC()*p1.getB()-p1.getC()*p2.getB(); 
			double zconstante=bdm/cbm;
			double yconstante=-p1.getC()/p1.getB()*zconstante-(p1.getD()/p1.getB()); 
			Pos3D lepoint=new Pos3D(0,yconstante,zconstante); 
			Pos3D ladirection=new Pos3D(1,0,0);
			return new DroiteSpace3D(lepoint,ladirection); 
		}
		
		double ab=p1.getA()*p2.getB()-p1.getB()*p2.getA(); 
		if(ab==0) // vecteur directeur dans le plan xOy (composante nulle selon z)
		{
			double bdm=p2.getB()*p1.getD()-p1.getB()*p2.getD(); 
			double cbm=p2.getC()*p1.getB()-p1.getC()*p2.getB(); 
			double zconstante=bdm/cbm;
			double yconstante=-1/p1.getB()*(p1.getC()*zconstante-p1.getD());
			Pos3D lepoint=new Pos3D(0,yconstante,zconstante);
			Pos3D ladirection=new Pos3D(1,-p1.getA()/p1.getB(),0); 
			return new DroiteSpace3D(lepoint,ladirection); 
		}
		double adm=p2.getA()*p1.getD()-p1.getA()*p2.getD(); 
		double acm=p2.getA()*p1.getC()-p1.getA()*p2.getC();
		double yconstante=adm/ab; 
		double zconstante=0; 
		double xconstante=-p1.getB()/p1.getA()*(adm/ab)-p1.getD()/p1.getA(); 
		double zcoef=1; 
		double ycoef=acm/ab; 
		double xcoef=-p1.getB()/p1.getA()*(acm/ab)-p1.getC()/p1.getA(); 
		Pos3D lepoint=new Pos3D(xconstante,yconstante,zconstante);
		Pos3D ladirection=new Pos3D(xcoef,ycoef,zcoef); 
		return new DroiteSpace3D(lepoint,ladirection); 
	}
	
	/** Calcule les angles de rotation amenant un plan parallele a un plan perpendiculaire a ox **/
	
	public CoupleDAngles rotations(Plane p){
		CoupleDAngles resu=new CoupleDAngles(); 
		if(p.getA()==0){
			if(p.getB()==0){
				resu.setAlpha(0); 
				if(p.getC()<0) resu.setBeta(Math.PI/2); else resu.setBeta(-Math.PI/2); 
				System.out.println("resu "+resu); 
				return resu;
			}
			else // A=0 et B!=0
			 {
			 if(p.getB()<0) resu.setAlpha(Math.PI/2); else resu.setAlpha(-Math.PI/2);
			 }
			 }
	    else // A!=0; 
	    	resu.setAlpha(-Math.atan2(p.getB(),p.getA())); 
			
		// On s'occupe de beta
		double leng=Math.sqrt(p.getA()*p.getA()+p.getB()*p.getB()); 
		resu.setBeta(-Math.atan2(p.getC(),leng));  
		System.out.println("resu "+resu); 
		return resu; 
		}
		
	
	
	
	public void run(String arg) {
		
		double mu=2; // r0 
		double c=1;  // r1
		double a=3;  // rc
		double b=Math.sqrt(a*a-c*c); 
		
		int nbCercles=10; 
		int nbDiv=10; // toujours pair ? 
		
		
		primeIP = new ColorProcessor(X,Y);
		
		// Tout tient dans une image carree
		double limit=2*(a+mu)*1.2; 
		double lmoins=-2*(a+mu)*1.2;
	
		setBounds(lmoins,limit,lmoins,limit); 
		setCoefs(); 
		
		
		
		Polygon background=new Polygon(new int[]{0,0,X,Y},new int[]{0,X,Y,0},4);
		primeIP.setColor(Color.WHITE);
		
		primeIP.fillPolygon(background);
		
		
		primeIP.setFont(new Font(null,Font.ITALIC,60));
		primeIP.setColor(Color.GREEN);
		primeIP.setLineWidth(3); 
		drawLine(-limit,0,limit,0); 
		drawLine(0,-limit,0,limit); 

		setProcessor("Verifying", primeIP);
		
		
		primeIP.setColor(Color.BLACK);
		primeIP.setLineWidth(5);
		
		
		
		// nbCercles/2 plans pour les cercles verticaux (theta)
		// nbDiv/2 plans pour les cercles presque horizontaux (alpha)
		// calcul des intersections. 
		Plane[] plansTheta,plansAlpha; 
		DroiteSpace3D[][] intersections; 
		
		plansTheta=new Plane[nbCercles/2]; 
		for(int i=0;i<nbCercles/2;i++){
			double th=2*i*Math.PI/nbCercles; 
			Pos3D p1=new Pos3D(calculX(th, 0),calculY(th,0),calculZ(th,0)); 
			Pos3D p2=new Pos3D(calculX(th, Math.PI/3),calculY(th,Math.PI/3),calculZ(th,Math.PI/3)); 
			Pos3D p3=new Pos3D(calculX(th, -Math.PI/3),calculY(th,-Math.PI/3),calculZ(th,-Math.PI/3)); 
			plansTheta[i]=new Plane(p1,p2,p3); 
		}
		
		plansAlpha=new Plane[nbDiv/2]; 
		for(int i=0;i<=nbDiv/4;i++){
			double alpha=2*i*Math.PI/nbDiv; 
			Pos3D p1=new Pos3D(calculX(0, alpha),calculY(0,alpha),calculZ(0,alpha)); 
			Pos3D p2=new Pos3D(calculX(Math.PI/3, alpha),calculY(Math.PI/3,alpha),calculZ(Math.PI/3,alpha)); 
			Pos3D p3=new Pos3D(calculX(-Math.PI/3, alpha),calculY(-Math.PI/3,alpha),calculZ(-Math.PI/3,alpha)); 
			plansAlpha[i]=new Plane(p1,p2,p3); 
		}
		for(int i=nbDiv/4;i<nbDiv/2;i++){
			double alpha=(i+nbDiv/4+1)*2*Math.PI/nbDiv;
			Pos3D p1=new Pos3D(calculX(0, alpha),calculY(0,alpha),calculZ(0,alpha)); 
			Pos3D p2=new Pos3D(calculX(Math.PI/3, alpha),calculY(Math.PI/3,alpha),calculZ(Math.PI/3,alpha)); 
			Pos3D p3=new Pos3D(calculX(-Math.PI/3, alpha),calculY(-Math.PI/3,alpha),calculZ(-Math.PI/3,alpha)); 
			plansAlpha[i]=new Plane(p1,p2,p3); 
		}
		
		intersections=new DroiteSpace3D[plansTheta.length][plansAlpha.length]; 
		for(int i=0;i<plansTheta.length;i++)
			for(int j=0;j<plansAlpha.length;j++)
				intersections[i][j]=inter(plansTheta[i],plansAlpha[j]); 
		
		// On va tracer les cercles et les coupures alpha....
		for(int i=0;i<nbCercles;i++){
			//dessiner un cercle avec ses coupures
			ImagePlus gamma=NewImage.createRGBImage ("cercle "+i,X,Y, 1, NewImage.FILL_WHITE);	
			ImageProcessor i1=gamma.getProcessor();
			FileSaver fs=new FileSaver(gamma); 
			i1.setColor(Color.BLACK);
			i1.setLineWidth(5);
			
			
			double theta=2*Math.PI*i/nbCercles; 
			double denom=a*a-c*c*Math.cos(theta)*Math.cos(theta); 
			double diamCarre=4*b*b*(c*Math.cos(theta)-mu)*(c*Math.cos(theta)-mu)/denom; 
			double rayon=0.5*Math.sqrt(diamCarre);
			
			double petitDenom=a-c*Math.cos(theta);
			
			double x0=(mu*(c-a*Math.cos(theta))+b*b*Math.cos(theta))/petitDenom; 
			double y0=b*Math.sin(theta)*(a-mu)/petitDenom; 
			
			double petitDenomP=a+c*Math.cos(theta);
			
			double xp0=(mu*(c+a*Math.cos(theta))+b*b*Math.cos(theta))/petitDenomP; 
			double yp0=b*Math.sin(theta)*(a+mu)/petitDenomP; 
			double dx=x0-xp0; 
			double dy=y0-yp0; 
			double verifRayon=0.5*Math.sqrt(dx*dx+dy*dy); 
			//drawLine(x0,y0,xp0,yp0,i1); 
			
			double vx=a*(mu*c*Math.sin(theta)*Math.sin(theta)+b*b*Math.cos(theta))/denom;
			double vy=b*Math.sin(theta)*(a*a-mu*c*Math.cos(theta))/denom;
			Vector2D center=new Vector2D(vx,vy);
			drawCircle(center,verifRayon,i1); 
			System.out.println(theta+" "+rayon+" **  "+verifRayon+" "+center); 
			// Les coupures
			Vector2D[] points=new Vector2D[nbDiv]; 
			for(int j=0;j<nbDiv;j++){
				double alpha=2*j*Math.PI/nbDiv;
				double denomi=a-c*Math.cos(alpha)*Math.cos(theta); 
				double xi=(mu*(c-a*Math.cos(alpha)*Math.cos(theta))+b*b*Math.cos(theta))/denomi;
			
				double yi=(b*Math.sin(theta)*(a-mu*Math.cos(alpha)))/denomi; 
				double zi=(b*Math.sin(alpha)*(c*Math.cos(theta)-mu))/denomi;
				Vector2D pointLocal=new Vector2D(xi,yi); 
				points[j]=pointLocal; 
				/*
				double localx=pointLocal.distance(center); 
				System.out.println("---->"+localx+"  "+zi); 
				if(xi>center.getX())localx=-localx;  
				points[j]=new Vector2D(localx,zi);
				*/ 	
			}
			for(int j=0;j<=nbDiv/4+1;j++){
				// tracer les lignes
				System.out.println(j); 
				drawLine(points[j],points[nbDiv/2-j],i1);
				if(j!=0)
				drawLine(points[nbDiv/2+j],points[nbDiv-j],i1);
			}
			
			// Ici, on devrait recalculer les coupures par intersections de plan (et translation pour remise en 0, et rotations...
			// waw le bordel
			/*
			for(int j=0;j<intersections[i].length;j++){
				// Calculer et afficher la jieme ligne de coupure
				DroiteSpace3D coupure=intersections[i][j];
				CoupleDAngles larot=rotations(plansTheta[i]); 
				DroiteSpace3D coupureRot=coupure.rotate(larot); 
				ArrayList<Sommet> lesBords=intersectCadre(lmoins,limit,lmoins,limit, coupureRot);
				if(lesBords.size()==2){ // tracer la ligne
					Sommet s1=lesBords.get(0); 
					Sommet s2=lesBords.get(1);
					drawLine(s1.getY(),s1.getZ(),s2.getY(),s2.getZ(),i1); 
				}// if lesBords
			}// for j
			
			*/
			
			gamma.show();
		}
		
		System.out.println("fini"); 
		
		
	
}
	
}
	